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ABSTRACT 

Self-consistent A-body simulations of modified gravity models are a key ingredient 
to obtain rigorous constraints on deviations from General Relativity using large-scale 
structure observations. This paper provides the first detailed comparison of the results 
of different A-body codes for the /(A), DGP, and Symmetron models, starting from 
the same initial conditions. We find that the fractional deviation of the matter power 
spectrum from ACDM agrees to better than 1% up to /c ~ 5 — lOhMpc”^ between the 
different codes. These codes are thus able to meet the stringent accuracy requirements 
of upcoming observational surveys. All codes are also in good agreement in their 
results for the velocity divergence power spectrum, halo abundances and halo profiles. 
We also test the quasi-static limit, which is employed in most modified gravity A- 
body codes, for the Symmetron model for which the most significant non-static effects 
among the models considered are expected. We conclude that this limit is a very good 
approximation for all of the observables considered here. 

Key words: Gosmology - (cosmology:) large-scale structure of Universe 


1 INTRODUCTION 

Cosmology, and in particular observations of the large-scale- 
structure (LSS), provide unique possibilities for testing Gen¬ 
eral Relativity (GR) on length scales that cannot be probed 
by any other means (see e.g. Jain & Khoury 2010; Koyama 
2015; Berti et al. 2015, for reviews). Motivated by the ob¬ 
served accelerating expansion of the universe (Riess et al. 
1998; Eisenstein et al. 2005; Bennett et al. 2013; Planck 
Collaboration et al. 2015), a number of theories have been 
proposed, in which the acceleration is explained by devia¬ 
tions from GR on large scales (Carroll et al. 2004; Dvali et al. 


2000; Nicolis et al. 2009; Hinterbichler & Khoury 2010), (see 
Clifton et al. 2012, for a comprehensive review). In almost 
all cases, these theories add to the standard massless spin-2 
graviton of GR a new light scalar degree of freedom if. Thus, 
the Einstein and fluid equations of standard cosmology are 
augmented by an equation of motion for (f. All these models 
are faced with the challenge of producing sizeable effects on 
large scales - the most desirable being a natural explana¬ 
tion of the accelerated expansion - while at the same time 
passing the stringent local constraints on modifications to 
GR (Will 2014). This requires some form of screening mech- 
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anism, that is, a way of dynamically suppressing the effects 
of the fifth force mediated by Lp in high-density regions (com¬ 
pared to the cosmological average) where local experiments 
have tested GR to high precision. Such an effect is typi¬ 
cally realised by nonlinearities in the equation of motion of 
(/?, which results in a violation of the superposition principle 
that suppresses the fifth force. The screening leads to a com¬ 
plex interplay between the large scale distribution of matter 
and the magnitude of the fifth force mediated by (p. 

The standard tools to compute large-scale structure ob¬ 
servables in the nonlinear regime are N-hodj simulations. 
Consequently, in order to robustly test GR with cosmology, 
reliable A/-body simulations of modified gravity models are a 
necessity. These simulations must solve the nonlinear equa¬ 
tion of p in conjunction with the Vlasov-Poisson system that 
is solved in standard iV-body simulations. The nonlinear na¬ 
ture of the scalar field equation requires the implementation 
of novel numerical techniques, which is what makes iV-body 
simulations of modified gravity so challenging. 

To date, a number of codes have been developed to per¬ 
form simulations of modified gravity. For instance, Oyaizu 
(2008) presented a code that simulates the Hu-Sawicki f{R) 
model (Hu V Sawicki 2007) by solving the scalar field equa¬ 
tion on a fixed mesh/grid (throughout we use the words 
mesh and grid interchangeably). Based on the work of Oy¬ 
aizu (2008), Schmidt (2009a,b) developed a code (which we 
call DGPM here) that performed simulations of the Dvali- 
Gabadadze-Porrati (DGP) braneworld model (Dvali et al. 
2000), also on a fixed mesh. Fixed grid simulations of the 
DGP model were also performed in Chan V Scoccimarro 
(2009) and Khoury V Wyman (2009). More recently, efforts 
have been made to simulate modified gravity cosmologies 
on adaptively refined meshes, which allow for better resolu¬ 
tion on small scales, where the effects of the screening are 
most important. These efforts resulted in the development 
of a modified version (Li V Zhao 2009; Zhao et al. 2011; Li 
et al. 2011; Li V Barrow 2011) of the mlapm code (Knebe 
et al. 2001), which is a serial iV-body code. The implemen¬ 
tation of modified gravity solvers on adaptive mesh refine¬ 
ment (AMR) parallelizable codes was achieved with the de¬ 
velopment of the ECOSMOG (Li et al. 2012, 2013a,b), mg- 
GADGET (Puchwein et al. 2013) and ISIS (Llinares V Mota 
2013; Llinares et al. 2014) A-body codes. While careful con¬ 
sistency checks have been performed by the authors of each 
code (for example, by solving test cases with known ana¬ 
lytical solutions), no detailed comparison between codes has 
been performed so far. A main goal of this paper is pre¬ 
cisely to provide a rigorous cross-check of the accuracy of 
the nontrivial algorithms of these codes. This is particularly 
important in light of the stringent accuracy requirements 
demanded by current and future observational campaigns. 

Here, we simulate the f{R), DGP and Symmetron (Hin- 
terbichler & Khoury 2010) models with the dgpm, egos- 
MOG, MG-GADGET and ISIS codes. We start all simulations 
of the different codes from the same initial conditions and 
compare their results for the matter and velocity divergence 
power spectra, halo mass function, as well as density, force 
and velocity profiles of dark matter haloes. For the case of 
the Symmetron model, we also measure the impact of assum¬ 
ing the quasi-static limit in A-body simulations of modified 
gravity, which amounts to neglecting time derivatives of p, 


by comparing with the results of a version of the ISIS code 
that explicitly solves for the time evolution of p. 

An important consideration in such a comparison 
project relates to determining the target accuracy. To guide 
ourselves in the interpretation of the results, we use the ex¬ 
pected accuracy of the next generation of LSS surveys. For 
instance, for a mission such as that to be carried out by the 
Euclid satellite^ (Laureijs et al. 2011; Amendola et al. 2013), 
the nonlinear matter power spectrum up to a wavenumber 
k ~ 5/iMpc“^ should be accurate to 1%^. Since our goal here 
is to accurately calibrate modified gravity effects, we aim for 
an agreement in the fractional change of the matter power 
spectrum relative to AGDM of 1% or better. The agreement 
of the different codes on the absolute AGDM predictions is 
not of primary concern here, as there are dedicated compar¬ 
ison projects for this purpose (Schneider et al. 2015). We 
shall also aim for an accuracy of a few per cent in the code 
results for velocity statistics, the halo mass function and 
halo profiles. 

Before proceeding, we note that the results presented in 
this paper contribute to recent efforts in testing and compar¬ 
ing A-body codes and codes that extract observables from 
simulations. For instance, there have been detailed com¬ 
parisons of standard GR A-body codes (Scannapieco et al. 
2012; Schneider et al. 2015) and also of codes that iden¬ 
tify dark matter haloes (Knebe et al. 2011), voids (Golberg 
et al. 2008), halo substructure (Onions et al. 2012, 2013; 
Pujol et al. 2014; Hoffmann et al. 2014), galaxies (Knebe 
et al. 2013a), tidal debris (Elahi et al. 2013), merger trees 
(Srisawat et al. 2013), halo mock generation (Ghuang et al. 
2014) and galaxy mass reconstruction (Old et al. 2014), just 
to mention a few. See Knebe et al. (2013b) for a review on 
the current status of structure finding in A-body simula¬ 
tions. Gomparison projects of these (often complex) numer¬ 
ical techniques are crucial to identify any worrying system- 
atics in the theoretical predictions. 

The rest of the paper is organised as follows. In Sec¬ 
tion 2, we briefly introduce and review the f{R), DGP and 
Symmetron models that are used in our numerical compar¬ 
ison. In Section 3, we outline the general numerical tech¬ 
niques that are used to solve the modified gravity equations 
in a A-body solver. In Section 4, we summarise the main 
features of each of the codes and specific details in how they 
tackle the equations of the different models. In Section 5, we 
present our comparison results for the matter power spec¬ 
trum (5.2), the velocity divergence spectrum (5.3), the halo 
mass function (5.4), and the density, force and velocity pro¬ 
files of haloes (5.5). We summarise our findings and draw 
our conclusions in Section 6. 


2 MODIFIED GRAVITY THEORY 

Alternative models to AGDM (like the modified gravity the¬ 
ories studied here) are numerous (see e.g. Amendola & Tsu- 
jikawa 2010; Glifton et al. 2012; Koyama 2015), as are the 
problems with which they must struggle. Some models are 

^ http://www.euclid-ec.org/ 

^ See for example Kitching & Taylor (2011) for requirements on 
the error envelope around the non-linear P{k) for cosmic shear 
tomography. 
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plagued by theoretical instabilities and others require at 
least some degree of fine-tuning of the model parameters 
in order to meet observational constraints. One particularly 
simple extension of GR is the inclusion of a single scalar 
field Lp to the standard GR Einstein-Hilbert action. However, 
when coupled to matter, the scalar field gives rise to an ad¬ 
ditional gravitational interaction, which is often referred to 
as a fifth force (see e.g. Amendola et al. 2013; Mota & Shaw 
2006; Hellwing & Juszkiewicz 2009; Hellwing et al. 2013a). 
This fifth force can be quantihed by 7 = |FFifth|/|FN| where 
Fn is the standard Newtonian gravitational force that we 
obtain in the weak-field limit of GR. Several experiments 
(see e.g. Adelberger 2002 ; Bertotti et al. 2003; Williams et al. 
2004; Will 2014) have constrained 7 <C 1 on Earth and in 
the Solar System. This seems to leave us with two possible 
explanations: either the fifth force is zero on all scales, i.e., 
7 = 0 , or 7 is not a constant but instead varies in space 
and/or time. 

Models where 7 is space-dependent are dubbed screened 
modified gravity models, as one typically desires the hfth 
force to be screened in high-density environments (like the 
Solar System). Next, we follow Joyce et al. (2015) in the clas- 
sihcation of the different types of screening mechanisms. A 
general Lagrangian density for the scalar field can be written 
schematically as 

c = 9^, - v{<p) + porp (i) 

where represents derivative self-interactions of the 

scalar field, V{lp) is a potential, is a coupling function 
and Tjf is the trace of the matter energy-momentum tensor. 
Eor nonrelativistic matter fields, Tjf = —Pm, the dynamics 
of if therefore depend on the local density of the system, pm- 
Around the background (p, the dynamics of the fluctuations 
of (p are determined by three parameters: the mass m{p) 
(roughly given by the curvature of the effective potential), 
the coupling and the kinetic function Screen¬ 

ing can be realised mainly in three different ways utilising 
these three parameters: 

• Large mass 

If the mass of the fluctuations is large in dense envi¬ 

ronments, then the scalar field cannot propagate beyond its 
Gompton wavelength and the fifth force mediated 

by the scalar field is suppressed. On the other hand, in low 
density environments such as the cosmological background, 
the mass can be light and the scalar field mediates a sizeable 
fifth force. This idea characterises the so-called Chameleon 
type of screening (Khoury & Weltman 2004a,b); 

• Large kinetic term 

If the kinetic function Z^^(p) is large in dense environments, 
the coupling to matter is suppressed. One can either make 
the hrst or the second derivative of the scalar field large 
in dense environments. The former case is realised in the 
k-mouflage (Babichev et al. 2009; Brax & Valageas 2014) 
and D-BIonic type of screening (Burrage & Khoury 2014), 
while the latter case characterises the Vainshtein screening 
mechanism (Vainshtein 1972); 

• Small coupling 

If the coupling to matter /3((p) is small in the region of high 
density, the strength of the fifth force Fpifth is weak and the 


modifications to gravity are suppressed. On the other hand, 
in low density environments, the size of the fifth force can 
be of the same order as standard gravity (7 ~ 1). This idea 
is realised in the dilaton (Brax et al. 2010) and Symmetron 
mechanisms (Hinterbichler V Khoury 2010). 

In this code comparison project, we take f{R), DGP 
and Symmetron gravity as our working example models that 
screen the fifth force via large mass, large kinetic terms and 
small coupling strengths, respectively. Thus, while we do not 
consider every individual modified gravity model proposed 
in the literature, our simulations do cover all classes of mod¬ 
els. Throughout, we work with the perturbed Eriedmann- 
Robertson-Walker (ERW) spacetime metric in the Newto¬ 
nian gauge 

ds^ = —(1 + 24>)dt^ + a^(l — 2'^)6ijdx^dx^, ( 2 ) 

where T and $ represent the two gravitational potentials. 
The dynamics of nonrelativistic matter is governed by 4>, 
whereas the bending of light is determined by the lensing 
potential 4>+ = (4> + T)/2. The modified gravity simula¬ 
tions employed here assume the same weak-field and non¬ 
relativistic limit as standard GR simulations, i.e. higher or¬ 
der terms in the dark matter velocities, as well as dynami¬ 
cally generated vector and tensor modes are neglected. 

In addition, unless otherwise specified, we assume the 
quasi-static limit for the modified gravity field equation. 
This refers to neglecting the time derivatives of the per¬ 
turbed helds, a.s in (p = (p 6p ^ (p, where is the fluc¬ 
tuation of the scalar field. In this paper, we shall assess the 
validity of the quasi-static limit in the A-body simulations 
of the Symmetron model. 

We now describe the specific modified gravity models 
considered in this paper. We will often refer back to the 
quasi-Newtonian potential which is defined through the 

Poisson equation, 

V^4>Ar = dTrGa^Spm , (3) 

where Spm is the matter density perturbation. 

2.1 f{R) gravity 

f(R) gravity is arguably the most well-studied modified 
gravity model in the nonlinear regime of cosmological struc¬ 
ture formation. In this model, one adds a function of the 
Ricci scalar R to the Einstein-Hilbert action 

S = + (4) 

where g is the determinant of the metric g^i, and Sm is the 
action of the matter fields fii. In the quasi-static and weak- 
held limits, the relevant equations for nonlinear structure 
formation can be written as 

^21 IGttG 2 c- 1 2 c-x-. /,-\ 

V 4> = —-—a 5pm + -a 5R^ (5) 

o u 

2 

V7fi = -^[SR + STrGSpm], (6) 

where Spm = pm — pm and 5R = R — R are the density 
and Ricci scalar perturbations, respectively (overbars denote 
background averaged quantities), and fn = df{R)/dR. In 
this formulation, fn plays the role of the scalar degree of 
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freedom if that determines the fifth force and that is solved 
for by the numerical codes. 

We specialise to the Hu-Sawicki model (Hu & Sawicki 
2007), which is characterised by 


f{R) 


fn 


2 Ci(i^/?Tl^)'^ 

TD -- 

( \ n — 1 

-R/m j 


( 7 ) 

( 8 ) 


where = H^Q^rn is a mass scale (not to be confused 
with the mass of the scalar fluctuations relevant for the 
Chameleon mechanism discussed above) and ci, C 2 and n 
are model parameters. Note that the field value is negative 
{fn < 0) which is necessary to ensure a positive mass of the 
scalar degree of freedom and hence stability of the theory. 
Since 


- ^ « 9.T^GPm - 2fiR) = 3m" (a”® + , (9) 

one recovers a ACDM expansion history by setting ci /c 2 = 
For values (Qrn,f^A) = (0.269,0.731) (as we con¬ 
sider in the simulations of this paper), then —R ^ and 
one can write 



( 10 ) 


For the simulations of this paper we always consider n = 1, 
and hence, the remaining free parameter is C 2 . However, in 
previous studies it has become more common to specify the 
Hu-Sawicki model not in terms of C 2 , but in terms of the 
equivalent value of Jr at the present day, Jrq. Note that by 
making use of Eqs. (9) and (10), one can eliminate 5R in 
favour of Jr in Eqs. (5) and (6). Eor completeness, we note 
that the modified Poisson equation, equation (5), can also 
be written as 

v"-l> = - ^V"/r. (11) 


This makes it explicit that in f(R) models the total gravita¬ 
tional force is governed by a modified gravitational potential 
-!> = -l>jv - I/r. 

We note that the modified gravitational equations de¬ 
fined above hold only for the dynamical potential of the 
model. The leasing potential in f{R) models (which are 
equivalent to scalar-tensor theories with a conformal cou¬ 
pling to matter) is not affected by the extra degree of free¬ 
dom (Brax et al. 2008) in the weak-field limit. 

The term SR on the right-hand side of equation (5) de¬ 
pends nonlinearly on Jr (cf. equation (10)). The nonlinearity 
is what gives rise to the Chameleon screening mechanism. 
The screening of the fifth force is determined by the depth 
of the gravitational potential $a 7 . A spherically symmetric 
object is screened if the thin shell condition 

IJroo - JrsI < ( 12 ) 

is satisfied, where fRoo is the Jr field away from the ob¬ 
ject and Jrs is that inside the object. In order to satisfy the 
Solar System constraint, the Milky Way galaxy with the po¬ 
tential \^n\ ~ 10“® needs to be screened. This imposes the 
constraint \fRo\ < 10if one assumes that the Milky Way 


galaxy is an isolated object in the cosmological background 
(Hu & Sawicki 2007). 

Eor the simulations presented in this paper, we consider 
models with |/i?o| = 10“^ (E5) and |/ho| = 10“® (F6). Al¬ 
though the former parameter value may already be in ten¬ 
sion with Solar System tests, we choose to simulate this 
model anyway, since it gives rise to larger fifth forces and 
places the screening threshold for halos at mass scales that 
are well resolved in the simulation^. Our main goal in this 
paper is to compare the different code predictions for the 
fifth force, and not so much to study the observational via¬ 
bility of the models. 


2.2 DGP 

The DGP model is an example of a braneworld model. In 
this model, matter is confined to live in a four-dimensional 
brane, embedded in a five-dimensional bulk spacetime. The 
action is given by 



where denotes the five-dimensional metric in the bulk, 
with R^^^ being the Ricci scalar for while g and R are 
the induced metric on the brane and its Ricci scalar, re¬ 
spectively. and G denote the five- and four-dimensional 
gravitational constants. The matter fields are confined 
to the four-dimensional brane. The relative sizes of the two 
gravitational strengths is a parameter of the model known 
as the crossover scale, rc. 


1 

“ 2 G 


(14) 


below which gravity looks four-dimensional, and above 
which the five-dimensional aspects become important. The 
cosmological solutions of this model are characterised by two 
branches of solutions. The normal branch requires a dark en¬ 
ergy term to be added to the four-dimensional part of the 
action to explain the accelerated expansion of the Universe 
(Sahni & Shtanov 2003; Lue & Starkman 2004; Schmidt 
2009b); the more appealing self-accelerating branch does not 
require a dark energy field, but it is in tension with CMB and 
supernovae data (Pang et al. 2008) and is also plagued by 
problems associated with the propagation of ghosts (degrees 
of freedom whose energy is unbounded from below) (Luty 
et al. 2003; Nicolis & Rattazzi 2004; Koyama 2007). In this 
paper, we focus on the normal branch of the DGP model. 
The dark energy component on the brane is adjusted to pre¬ 
cisely yield a flat AGDM background cosmology (Schmidt 
2009b). 

The modifications to the gravitational law in this model 
are determined by a scalar field, (^, which is associated with 
the bending modes of the 4D brane. The brane-bending 
mode influences the dynamics of particles through the dy¬ 
namical potential 4>, which, assuming the same boundary 


^ See Fig. 6 in Gronke et al. (2015b). For |/ho| = 10“^ halos 
with mass M > 3 x 10^^ Mq//i are screened while the smallest 
halos we can resolve have M ~ 10^^ Mq//i. 
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conditions for $ and is given by 
$ = + 2^' 

The equation for ip reads, in the quasi-static and weak-field 
limits (Koyama & Silva 2007), 

(16) 

where (ViVj(y^)^ = (V^ Vj(/?)(V*V-^V), and the function (3 {a) 
is given by 


/?(a) = l + 2if(a)re(^l + ^|^) , (17) 

where we have assumed the normal branch of the DGP 
model already. The quasi-static approximation in the DGP 
models was tested for self-consistency in Section IV G of 
Schmidt (2009a) and was also recently shown to be an excel¬ 
lent approximation in Brito et ah (2014); Winther V Ferreira 
(2015). 

As in f{R) gravity, the propagation of photons, deter¬ 
mined by the lensing potential 4>+, is not directly affected by 
ip. In models like DGP it is the Vainshtein screening mech¬ 
anism that provides the chance to pass Solar System tests. 
For simplicity, we focus on spherically symmetric configu¬ 
rations to illustrate how the screening works. Writing down 
equation (16) in spherical coordinates and integrating once 
as J r^dr, one gets 


5 ^ I f GM{r ) 

3/3 v r / V r / 3/3 


(18) 


where M(r) is the mass enclosed inside a radius r, and a 
comma denotes partial differentiation. The solution to the 
last equation is given by 





GM{r) 


where we define the distance scale 


rv(r) 


n&rlGM{r)V'^ 

\ ) 


(19) 


( 20 ) 


which is known as the Vainshtein radius. This radius defines 
the distance from the centre of the spherical overdensity be¬ 
low which the spatial gradient of p becomes suppressed (and 
hence the fifth force effects become negligible). Explicitly, 
for a top-hat density profile of radius Rth and mass Mth, if 
r rv > Rth, then 


p,r _ 1 GMth _ 1 

~ ^ r 2 “ ^ 


( 21 ) 


i.e., the fifth force becomes a sizeable fraction of the force in 
GR (cf. equation (15)). On the other hand, if Rth < r rv, 
then p,r 0. 

In the simulations of this paper, we consider two param¬ 
eter values, PcHo = 1 and TcHo = 5. These were chosen to 
roughly match the F5 and F6 models, respectively, in terms 
of the values of erg at z = 0. 


2.3 Symmetron 


Olive V Pospelov (2008)), whose action is given by 


S' = 


+ 



Smi^gpiv, V^)- 


R 

IGttG 


-VO 


( 22 ) 


The matter fields, -0, couple to the Jordan frame metric 
which is given by a conformal rescaling of the Einstein frame 
metric g^iy 

gpbv — A {p)gfxu- (23) 


In the Symmetron model, the coupling function A{p) is 
given by 


where M is a mass scale. This coupling function determines 
the fifth force and the total gravitational force is given by 


F + 2 


— T 


pVp 
M 2 ■ 


(25) 


The potential is taken to be of the symmetry breaking form 

VO = H + jA/. (26) 

With these choices for A{p) and V{p)^ the model becomes 
invariant under the symmetry p ^ —p. The value of Vo is 
determined by the condition that the model gives rise to the 
observed accelerated expansion of the Universe (Hinterbich- 
ler et al. 2011). The field equation for p follows from the 
variation of the action, equation (22), with respect to p and 
reads 


\Z\p — 


(27) 


where, for nonrelativistic matter, the effective potential is 
given by 

VOv) = Vo+\[^-p)P + \\P. (28) 


In working with the model, it is convenient to define a matter 
density scale for symmetry breaking, pssB, and its associated 
scale factor, assB, where p(assB) = pssB, as 


U1 

U1 

td 

III 

to 

II 

CO 

b 

na/dsSB- 

(29) 

Other useful quantities are 



PoMy>\ . 1 

" M2 ’ " 72a* ’ 

P 

(30) 


where fio is the coupling strength of the unscreened fifth 
force, Ao is the Gompton length (giving the range of the fifth 
force) and po is the symmetry breaking vacuum expectation 
value (VEV) of p (when pm = 0). In the simulations of this 
paper we consider /3o = 1, Ao = l/i“^Mpc and assB = 0.5 
which lie on the boundary of the allowed {assB, Ao} param¬ 
eter space coming from local constraints (see the discussion 
below equation (32)). These parameter values have been pre¬ 
viously simulated in Davis et al. (2012); Brax et al. (2012a); 
Llinares et al. (2014); Llinares & Mota (2014). 

In the quasi-static limit, equation (27) becomes 


T—72 

V X 


f Pm 
2 Ao VPssB 



(31) 


The third model that we consider is the Symmetron model where x = The full field equation, without applying 

(Hinterbichler & Khoury 2010) (see also Pietroni (2005); the quasi-static limit, is discussed in Section 4.2.3. 
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During the cosmological evolution (Brax et ah 2011b, 
2012b; Hinterbichler et ah 2011; Davis et ah 2012) the field 
sits close to the global minimum of the effective potential 
at (/? = 0 for a < assn- For a > assB the effective potential 
develops two minima at (/? = ±(fo^yl — clqqb/ci^, to which 
the field at (/? = 0 (now a maximum) evolves, thereby spon- 
tanously breaking the (f ^ —(p symmetry. Since the field can 
choose different minima (+ or — branches) in different parts 
of the Universe, this model therefore leads to the forma¬ 
tion of domain walls. The properties of these domain walls 
have been studied beyond the quasi-static approximation 
(Llinares & Mota 2013; Llinares & Pogosian 2014; Pearson 
2014). 

Screening in the Symmetron model is very similar to 
that in the Chameleon//(i?) cases in the sense that the con¬ 
dition for screening is determined by the local gravitational 
potential. There is however the important difference that 
the coupling /3{(p) = which is constant for f{R) grav¬ 
ity, now depends on the local field value. In high density 
regions, pm > Pssb, the field moves towards (/? = 0, and 
since the coupling is proportional to the fifth force is 
suppressed. In the Symmetron model, the condition for the 
thin-shell effect is given by 

< ^AT/JCv’oo), (32) 

^0 

where poo is the p field far away from the object and ps 
is that inside the object. In order to satisfy Solar System 

bounds, we get the constraint (^ /^-iMpc ) ^ssb ~ ^iX) t>y 
assuming that the Milky Way galaxy is an isolated object 
in the cosmological background (Brax et al. 2012b). 

3 MODIFIED GRAVITY SIMULATIONS 
3.1 General force calculation 

Cosmological dark-matter iV-body simulations for standard 
gravity are characterised by the following two equations. 
First, we have the Poisson equation (3) 

= A-KGaSprr^, (33) 

which determines the Newtonian potential, 4>Ar, given the 
density fiuctuations 5pm - Second, we have the geodesic equa¬ 
tion 

x + 2i7x = -V4>, (34) 

which tells the particles how to move. At every timestep in 
the simulation, one (i) computes the density field from the 
particle positions; (ii) uses it in the Poisson equation to solve 
for the potential; and (hi) plugs 4> into the geodesic equation 
to move the particles. This process is repeated from some 
initial redshift, z — zi^ until typically z = 0. 

As we have seen in the previous section, modified grav¬ 
ity models alter this picture by modifying the Poisson equa¬ 
tion that governs the total gravitational potential (like in 
f{R) and DGP gravity), or by adding extra terms to the 
right-hand side of the geodesic equation (like the term 
(X VA{ip) in the Symmetron model)^. In general, modified 
gravity models can also differ in the background expansion 

Note that any extra term in the geodesic equation can always be 


rate of the Universe. In this paper, however, we shall always 
consider models with the same ACDM background evolu¬ 
tion. 

The bulk of the computing time in modified gravity sim¬ 
ulations is spent solving the nonlinear equation that governs 
p (see the previous section), which can be cast in the general 
form 

L[p\ = S{5pm,p), (35) 

where L is some nonlinear operator that acts on c^, and S 
is a source term that depends on the matter density fiuctu¬ 
ations and possibly on the scalar field. The exact functional 
form of L and S varies from theory to theory, but as we 
have discussed in the previous section, this operator should 
possess some degree of nonlinearity to ensure the presence 
of screening effects. The nonlinearity in the equations, how¬ 
ever, is what makes A-body simulations of these models so 
challenging. On the other hand, equation (33) is a linear el¬ 
liptic partial differential equation (PDF), which means that 
it can be solved with efficient fast Fourier transform (FFT) 
methods. This is in general not possible in modified grav¬ 
ity models, which typically have nonlinear equations. This 
difficulty can be overcome by employing a FFT-relaxation 
method (Chan & Scoccimarro 2009), if the equations are 
to be solved on a regular grid. However, this method does 
not work on irregularly-shaped refinements. The codes we 
compare in this study solve equation (35) via direct dis¬ 
cretisation and relaxation on such an irregular grid. 

In the rest of this section, we briefiy outline the relax¬ 
ation algorithm, describing also the main idea behind multi¬ 
grid acceleration methods. The latter significantly improves 
the efficiency of the relaxation algorithms. 

3.2 Iterative methods with multigrid acceleration 

With the exception of the isis-nonstatic code (see next sec¬ 
tion), all codes make use of multigrid acceleration to speed 
up numerical convergence of the partial differential equa¬ 
tions for p. Here, we briefiy review the main aspects of these 
techniques and refer the reader to Brandt (1977); Wesseling 
(1992); Trottenberg et al. (2000) and to the code papers 
(cf. Table 1) for further details on their implementation. 

3.2.1 Gauss-Seidel iterations 

The goal is to solve a differential equation that can be writ¬ 
ten in the form of equation (35). The basic algorithm con¬ 
sists in discretising the equation on a grid and using an 
iterative scheme to obtain improved solutions given an ini¬ 
tial guess. All codes assume periodic boundary conditions 
on the domain (unrefined) grid. The codes that include grid 
refinements use fixed boundary conditions on the boundary 
of the refinements, obtained by interpolating from the next 
coarser refinement level. Upon discretisation, the solution 
to equation (35) is given by the solution of the large set of 
algebraic equations 

L\P] = S\ (36) 

absorbed into the definition of a modified gravitational potential. 
This essentially illustrates the equivalence between the Jordan 
and Einstein frames. 
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where P and are the discretised versions of the L and S 
operators and is the field solution we aim to determine. 
The index I labels the refinement level of the grid. The dis¬ 
cretisation of the equation consists in writing each of the 
derivatives that appear in L as a combination of the values 
of (f on the grid cells. For instance, the codes employed here 
use the 3- and 4-point stencils 
2 1 

dxdyifij^k — + + — r^i-\-l,j-l,k 

^i — l,j-\-l,k ^i — l,j — l,k) 5 ( 38 ) 

where labels each grid cell. The iterations can be 

made in two different ways. If the operator L is linear, then it 
is best to perform explicit iterations, in which one rearranges 
the discretised equation analytically to solve directly for ipijk 
in each cell (this is, for instance, how the standard ramses 
code solves the Poisson equation). For nonlinear problems 
(as those in modified gravity), an implicit iteration scheme is 
more suitable. In this case, one can use the Newton-Raphson 
method to solve the equation 

T'[v5*] = - S'* = 0. (39) 


The resulting value of (p can be written as 

rpl 


(40) 


where the barred field corresponds to the updated value. Fi¬ 
nally, one needs to specify the way in which sweeps are made 
across the grid in each iteration step. The simplest sweep¬ 
ing strategy consists in following a lexicographic ordering, 
in which the calculation on the cell {z,j,/c} is followed by 
the calculation on the cell {i + 1, j, /c}, and so on. A more 
widespread strategy which has better convergence and par¬ 
allelisation properties is the so-called two colours scheme, in 
which the calculation is performed alternatively in cells of 
the same colour^ as in the colour scheme of a chess board. In 
the first half-sweep all black cells are updated, and the sec¬ 
ond half-sweep takes care of the remaining white cells. Fur¬ 
ther generalisations of this scheme exist with four or even 
eight colours. 

The iterations proceed until a certain convergence cri¬ 
terion is fulfilled. There are several criteria in the literature. 
These typically involve computing the residual, e\ of the 
solution defined as 


e* =L* [¥?*]-5'*. (41) 

The convergence criterion is then given by || e* ||< Cconverged, 
where || . || is a norm (typically L 2 ) that is taken over the en¬ 
tire grid and Cconverged is a (small) predefined constant that 
we use to define convergence. Exact solutions of the algebraic 
equation have = 0. However, owing to truncation errors 
(the error inevitably introduced by discretising a continuous 
equation), exact solutions of the algebraic equation are not 
equal to the solutions of the discretised equation. The itera¬ 
tions are then assumed to have converged once the residual 
falls below a predetermined fraction of the truncation error. 


However, the convergence becomes considerably slower as 
one approaches the true solution. This slowdown of the con¬ 
vergence is attributed to the fact that the components of the 
residual whose Fourier wavelength modes are longer than 
the size of the grid cell decay much more slowly than those 
modes whose wavelength is comparable to the grid size. The 
goal of multigrid methods is to speed up the convergence 
of these longer wavelength modes by using a hierarchy of 
coarser grids. In short, when the convergence of the solution 
starts to slow down, one interpolates the equation onto the 
next coarser grid of the hierarchy and solves it there. This 
makes the longer wavelength modes decay faster, therefore 
bringing the solution closer to its true value. This coars¬ 
ening scheme can proceed up to several coarser grids. The 
coarser solutions can then be interpolated back to the finer 
(original) level. 

To be more concrete, the typical way to arrange differ¬ 
ent resolutions is to use a set of grids whose resolutions are 
half, one quarter and so on of the target resolution. A two 
grid scheme is defined in the following way. One starts by 
performing a given number of iterations on the target grid 
1. Then, when the convergence slows down, one moves to 
the next coarser grid I — 1 and performs iterations for an 
equation whose solution corresponds to the error of 

the previous solution. In the case of a linear PDE (we will 
turn to the nonlinear case below), the equation that must 
be solved on the coarser grid is 

C-i[5^*-i]=i?(e*), (42) 

where L is now a linear operator while R is a restriction 
operator that is chosen according to the problem and trans¬ 
lates information from the fine grid (the target grid) to the 
coarse grid. Once the coarse grid iterations are done, one 
corrects the fine grid solution as 

1 ^* = (p* + T’((5(/3*“^), (43) 

where P is now a prolongation operator which translates 
information from the coarse to the fine grid. In general, one 
uses more than one coarse grid and these processes of going 
up and down in resolution are called V-cycles. After one V- 
cycle, if convergence is not yet achieved on the target grid, 
then further V-cycles are performed. All the codes analysed 
here use V-cycles, although other arrangements are possible 
such as W-cycles (in these, one can move in between coarser 
levels several times before returning to the target grid). 

In the case of nonlinear equations however, this multi¬ 
grid algorithm requires some changes, as it relies on the lin¬ 
ear superposition of solutions from different grids. In the 
nonlinear case, instead of solving following the solution for 
the errors on the coarse grids, one obtains improved approx¬ 
imations of the solution (not the error) itself. In this case, 
the coarse grid iterations are made according to 

L'- 1 [VP*-'] = -R{e‘ (v*, 5*)) + e*-' {Rp), Rp)) (44) 

and the coarse grid correction of the fine grid solution is 
given by 

<p'-= - Rp)). (45) 


3.2.2 Multigrid acceleration 

In a relaxation method such as the one described above, dur¬ 
ing the first few iterations the residual decays very efficiently. 


4 CODE AND ALGORITHM DESCRIPTION 

In this section, we briefly introduce the different A-body 
codes compared in this paper and comment on some aspects 
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of the numerical handling of the specific model equations. 
We shall keep our description simple and refer the interested 
reader to the code papers for the details. Some of the key 
features of the codes are summarised in Table 1. 

4.1 Code summary 

4 . 1.1 DGPM 

The DGPM code (Schmidt 2009a, b) is a fixed-grid particle- 
mesh code that solves the Vainshtein-type equation of mo¬ 
tion (16). Based on the fixed-grid f{R) code presented in 
Oyaizu (2008), it employs a second-order leapfrog scheme 
with fixed step size Aa in scale factor to advance parti¬ 
cles. Densities are interpolated onto the grid using cloud- 
in-cell (CIC) interpolation, which is also used to evaluate 
derivatives on the grid. The Poisson equation for the New¬ 
tonian potential is solved using FFT on the fixed grid. The 
Gauss-Seidel relaxation is then performed using the Newton- 
Raphson method and multigrid acceleration as described in 
Section 3.2. At each multigrid level, 5 to 10 relaxation sweeps 
are performed. Gonvergence to an RMS residual of less than 
10 “^° (where typical values of (p are of order 10“^) is usually 
reached within 3 V-cycles. The step size used for this paper 
is Aa = 0.02, which results in 490 steps from z — 49 (when 
all our simulations start) to = 0. 


4 . 1.2 ECOSMOG 

The EGOSMOG code (Li et al. 2012) is built on top of the 
publicly available adaptive mesh refinement (AMR) A-body 
code RAMSES (Teyssier 2002). The code can be compiled 
to work with GIG (like in ramses) or triangular shaped 
cloud (TSG) schemes for the interpolation of the density 
and force fields. Unless otherwise specified, the egosmog 
results shown in this paper are for TSG. The time evolution 
is performed with a second-order leapfrog algorithm with 
adaptive timesteps (set by the AMR grid, as in ramses). 
The Gauss-Seidel relaxations of the scalar field equation are 
performed using the Newton-Raphson methods with multi¬ 
grid acceleration on all levels of the AMR grid. For the simu¬ 
lations in this study, the grid was set to be refined whenever 
the particle number exceeded 8 inside a grid cell. 

In addition to the f{R) and DGP gravity results shown 
in this paper, the egosmog code has also been used to sim¬ 
ulate dilaton (Brax et al. 2011a), Symmetron (Davis et al. 
2012 ; Brax et al. 2013), Gubic Galileon (Barreira et al. 
2013a) and Quartic Galileon (Li et al. 2013b) gravity cos¬ 
mologies, as well as a general parametrization of Ghameleon 
theories (Brax et al. 2012a). Different versions of the code 
differ in the detailed way the scalar field equations are 
solved. The performance of the relaxation algorithm is also 
slightly different, although, for all these models, the resid¬ 
uals always reach a value of < 10“^ times the truncation 
error after 5—10 V-cycles. The extensions made to ramses 
to develop egosmog can also be straightforwardly coupled 
to the hydrodynamic modules of the base code, although to 
date such a project has never been undertaken. In case of 
the DGP model, egosmog solves a different version of the 
scalar field equation as dgpm, which will be discussed in 
Section 4.2.2 below. 


4 . 1.3 MG-GADGET 

The MG-GADGET code (Puchwein et al. 2013) is an exten¬ 
sion of the cosmological hydrodynamical TreePM+SPH sim¬ 
ulation code P-GADGET3 which is itself based on gadget2 
(Springel 2005). It features the baryonic physics modules of 
P-GADGET3 as well as a modified gravity solver. In addi¬ 
tion, MG-GADGET allows the inclusion of massive neutrinos 
in simulations of modified gravity (Baldi et al. 2014) by mak¬ 
ing use of the particle-based massive neutrino module (Viel 
et al. 2010) which is implemented in p-gadget3. 

In contrast to ramses, p-gadget3 does not intrinsi¬ 
cally possess an AMR grid. To overcome this, mg-gadget 
constructs an adaptively-refining grid that covers the whole 
simulation volume by appropriately choosing nodes from the 
oct-tree structure of p-gadget3’s Poisson solver. This grid 
is then used to solve for the scalar degree of freedom using 
the method described in Section 3.2, i.e. using GIG den¬ 
sity assignment and multigrid-accelerated Gauss-Seidel re¬ 
laxation on the different levels of the AMR grid. 

So far the code has been used to simulate the Hu V 
Sawicki (2007) f{R) gravity model, both in collisionless (DM 
only) and hydrodynamical simulations (Arnold et al. 2014, 
2015). For hydrodynamical simulations the fluid equations 
are solved using the same entropy-conserving SPH scheme 
(Springel V Hernquist 2002) as p-gadget3. 


4 . 1.4 ISIS 

The ISIS code (Llinares et al. 2014), like egosmog, is a mod¬ 
ified version of ramses. To date, isis has been used to sim¬ 
ulate f{R) gravity (Llinares et al. 2014), the Symmetron 
model in both the quasi-static (Llinares et al. 2014) and 
non-static limits (Llinares V Mota 2014), the non-static dis- 
formal gravity model (Koivisto et al. 2012) in its pure disfor- 
mal limit (Llinares V Mota 2015), the non-static disformally 
coupled Symmetron model (Hagala et al. 2015) and the Gu¬ 
bic Galileon / DGP model (Winther V Ferreira 2014). In 
Hammami et al. (2015) and Hammami V Mota (2015), ISIS 
has also been used to study hydrodynamic effects in simu¬ 
lations of f{R) and Symmetron models. 

The static version of ISIS solves the equation of motion 
of the scalar field using the multigrid methods outlined in 
the previous section. The code uses a GIG scheme to inter¬ 
polate the density from the particles to the grid, and the 
time steps of each particle are determined by the AMR grid 
(as in standard ramses). In our simulations of the static 
ISIS code, each grid cell was refined whenever the particle 
number contained in it exceeded 8 (as in the egosmog sim¬ 
ulations) . 

The non-static version of ISIS (the version that goes be¬ 
yond the quasi-static limit) uses a leapfrog scheme to evolve 
the scalar field in time. In this case, the code includes two 
time steps: a coarse one for the particles, which is deter¬ 
mined by the domain grid (the non-static version does not 
admit refinements) and a finer one for the time evolution of 
the scalar field. The total number of scalar field time steps 
is typically three to four orders of magnitude larger than the 
number of particle time steps. 
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Table 1. Key features of the A/-body codes compared in this paper 


Code 

DGPM 

EGOSMOG 

MG-GADGET 

ISIS 

ISIS-NONSTATIG 

Code paper 

Schmidt (2009a) 

Li et al. (2012, 2013a) 

Puchwein et al. (2013) 

Llinares et al. (2014) 

Llinares Mota (2014) 

Base code 

Oyaizu (2008) 

RAMSES 

P-GADGET3 

RAMSES 

RAMSES 

Density assignment 

GIG 

GIC/TSG 

GIG 

GIG 

GIG 

Force assignment 

GIG 

GIC/TSG 

effective mass 

GIG 

GIG 

Adaptive refinement? 

No 

Yes 

Yes 

Yes 

No 

Timestep 

Fixed 

Adaptive 

Adaptive 

Adaptive 

Adaptive 

MG solver 

Multigrid 

Multigrid 

Multigrid 

Multigrid 

Leapfrog 

Gravity solver 

Multigrid 

Multigrid 

TreePM 

Multigrid 

Multigird 

Parallelisation 

OpenMP 

MPI 

MPI 

MPI 

MPI 

Programming language 

C + + 

Fortran 

C 

Fortran 

Fortran 

Models simulated 

DGP 

/(R)/DGP 

RR) 

/(i?)/Symmetron 

Symmetron 


4.2 Model algorithms 

In the remainder of this section, we briefly outline the strat¬ 
egy employed by the different codes to solve the equations 
of the three modified gravity models we consider. 


4-2.1 f{R) simulations 

The simulations of the f{R) Hu-Sawicki model were per¬ 
formed with ECOSMOG, MG-GADGET and ISIS. All these codes 
discretise and relax the scalar field equation of motion, equa¬ 
tion (6), in a similar way. Instead of solving for fn directly, 
the scalar field is redefined in terms of u = ln(/i?//i?(a)) 
which is then numerically computed. Using the variable u 
has considerable advantages in terms of numerical stability 
as it implicitly avoids unphysical positive values of fn when 
performing the Newton-Raphson iterations. 

Once fn is found, the three codes compute the total 
force in slightly different ways. In egosmog the code uses 
the solution for fn to compute the SR term on the right- 
hand side of equation (5). The code then determines the 
total potential 4> by solving the modihed Poisson equation 
in a similar way to the standard gravity solver in ramses. 
The total (modihed) force, V4>, is hnally interpolated from 
the mesh to the particle positions (like in standard ramses). 

In mg-gadget, the total force is also obtained by solv¬ 
ing the modihed Poisson equation, but by making use of 
the standard tree+particle-mesh gravity algorithm of the 
base code. In order to do so, the modihed Poisson equa¬ 
tion is rewritten in terms of an effective mass density: 
V4> = 4:7TG(Sp-\-Spef¥), where (5peff = 24 ^^^- Adding 

the effective to the real mass density, the values of the scalar 
held can be directly used in the highly optimized and effi¬ 
cient TreePM gravity algorithm of p-gadget3 to compute 
the total force. 

Although these two methods for solving the modihed 
Poisson equation are mathematically equivalent, they can 
yield different numerical accuracies. In addition to the con¬ 
venience of using the standard Poisson solver, mg-gadget’s 
effective mass algorithm avoids an interpolation of the scalar 
held gradient (the hfth force) from the adaptive mesh to 
the particle positions, as the tree force is directly computed 
there. It might, nevertheless, be somewhat less accurate in 
highly screened regions than the method used in egosmog 
due to numerical summation errors in the tree gravity, which 



Figure 1. Fifth to Newtonian force ratio profile, ob¬ 

tained from interpolating the scalar field gradient from the grid 
to the particle positions (blue) and from the effective matter den¬ 
sity method using the tree force calculation (red) in mg-gadget 
for the f(R) simulations. The profiles shown are the average for 
haloes with masses M G [1 x 5 X Mq/Zi. The errorbars 
are the variance of this average. 

can result in less precise screening of the hfth force. This 
causes somewhat larger random force errors for the indi¬ 
vidual particles while the integrated effects are expected to 
average out. Fig. 1 displays prohles of the ratio of hfth to 
Newtonian force in dark matter haloes obtained via inter¬ 
polation of the gradient from the grid (blue) and via the ef¬ 
fective density scheme (red). Although it is noticeable that 
there are signihcant differences between the two methods 
on small scales, this occurs only in a regime where the hfth 
force is already highly screened (< 1% of normal gravity). 
As force errors of around one percent also occur in the stan¬ 
dard tree gravity algorithm, these errors are expected to be 
negligible. In the region where screening just sets in (which 
several observables might be sensitive to) the curves almost 
perfectly match each other. Consequently, the mentioned in¬ 
accuracies in the hfth force calculation will not change the 
total force signihcantly and will therefore only have a very 
minor impact on observables (as we shah see in more detail 
in the next sections). 
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Finally, in isis, the code computes the gradient of the 
scalar field and interpolates it onto the particle positions. 
The GR gravitational potential, is solved exactly as in 
RAMSES, and the total force at the particle positions is given 
by — V/i ?/2 (cf. equation ( 11 )). 

We refer the reader to the respective code papers for 
more details about how ecosmog, mg-gadget and isis 
solve the modified gravity equations, code tests, and explicit 
discretisation of the equations. 


4 . 2.2 DGP simulations 


The simulations of the DGP braneworld model shown here 
were performed with the dgpm and egosmog codes. In 
DGPM, the scalar field equation is discretised as it is written 
in equation (16). Some of the simulations using DGPM pre¬ 
sented in Schmidt (2009a,b) employed a Gaussian smooth¬ 
ing of the density field in order to improve convergence. For 
the analysis of this paper, however, no such smoothing was 
performed. 

In the strategy employed by the egosmog code, equa¬ 
tion (16) is manipulated analytically before being discre¬ 
tised. This manipulation is called the operator-splitting trick 
(Ghan & Scoccimarro 2009; Li et al. 2013a,b), which we de¬ 
scribe next. Equation (16) can be cast as 


(1 — w) (V^(/?)^ + aV‘^ip — E = 0 , 


where 


3/3(a)a^ 


(46) 

(47) 


E = + (48) 

and w is a constant numerical factor. Equation (46) can be 
solved once to yield 


2 a ± + 4(1 - w)S 

- 2 ( 1 ^) -• 


(49) 


Decomposing the term ViVj(p into its trace and traceless 
part (this is the operator-splitting trick), 

ViVjtp = + ViVj(p, (50) 

it is possible to show, after a bit of algebra, that 

E=(VzVj(/?) H-^— Sp, if w = (51) 

That is, if w = 1/3, then this cancels out the term in 
E. The reason why this is useful is because, upon discreti¬ 
sation, the traceless part ViVjp does not depend on the 
grid cell value, pijk, but only on its neighbours. Hence, in 
the equation that egosmog solves, which is equation (49), 
(fijk appears only on the left-hand side, and not inside the 
square-root. This is found to improve considerably the per¬ 
formance of the code. Moreover, by taking out pijk from 
inside the square-root one also avoids potential problems 
associated with imaginary square-roots caused by some bad 
initial guess for pijk. The sign of the square-root in equa¬ 
tion (49) is chosen to be the same as the sign of the a func¬ 
tion. This is the solution which corresponds to the physical 
(linear theory) result that ^ 0 , when ^ 0 . 

Once (f is found on every grid cell, both dgpm and 


egosmog compute the total force, V4>, as the sum of normal 
gravity and the fifth force, V4> = V4>a( + V(/?/2. 

Eor completeness, we point out that if Sp becomes neg¬ 
ative (as it does in voids), then there is the risk that the 
argument of the square-root in equation (49) may become 
negative. This does not happen for the DGP model, but sim¬ 
ilar Vainshtein screening models such as the Gubic (Barreira 
et al. 2013a) and Quartic Galileon (Li et al. 2013b; Barreira 
et al. 2013b) do suffer from imaginary square-root problems 
in low-density regions (see Winther & Eerreira (2015) for a 
discussion about the meaning of these imaginary square-root 
problems). 


4 . 2.3 Quasi-static and non-static simulations of the 
Symmetron model 


The simulations of the f{R) and DGP models are performed 
under the quasi-static approximation. To go beyond this ap¬ 
proximation means to explicitly take into account the time 
derivative terms of the scalar field in the equations. This 
way, the solution for the scalar field at a given time depends 
also on its past evolution, as opposed to depending only on 
the matter configuration at that given time. To date, non¬ 
static cosmological simulations of modified gravity have been 
performed for the Symmetron and disformal gravity models 
using the explicit leap-frog method (Llinares V Mota 2013, 
2014; Hagala et al. 2015), for f{R) gravity using the im¬ 
plicit Newton-Gauss-Seidel method (Bose et al. 2015) and 
for the DGP and the Gubic Galileon models using both of 
the methods mentioned above, but only in a spherical sym¬ 
metric spacetime (Winther & Eerreira 2015). 

Eor the Symmetron model, the full Klein-Gordon equa¬ 
tion (27) reads 


X + 3iLx- 



1 


%x^-x + x^' 


(52) 


where 77 is the matter density in units of the background 
value. In isis-nonstatig (the modified version of isis that 
relaxes the quasi-static approximation), the second order 
equation of motion of the scalar field is decomposed into 
a system of two first order equations as 


X 

Q 


a'^ 


V72 ^ 


%x^-x + x^' 


(53) 

(54) 


which are used to propagate both x and q using a leapfrog 
algorithm. The "position" x and "velocity" q are displaced 
from each other by 1/2 timestep and the discretised equa¬ 
tions become 


Xn = Xn-1 +Xn-l/2At, (55) 

qn+ 1/2 = qn- 1/2 + At, (56) 

where fn = f{tn). The spatial derivatives of x in the formu¬ 
las above are calculated from the grid using a 5-point stencil, 
see Llinares V Mota (2014) for a detailed description of the 
implementation of the scheme. 

The quasi-static simulations of the Symmetron model 
are performed by discretising equation (52) and neglecting 
the first two terms on the left-hand side. In this case, one 
does not need to explicitly evolve the scalar field. Instead, 
given the matter distribution, 77 , at a given time step, the 
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code relaxes the equation that contains only the scalar field 
and its second spatial derivative (not time derivatives), using 
the multigrid methods described in Section 3.2. 


5 RESULTS 

In this section we present the main results of this code com¬ 
parison project. We present and discuss the different code 
results for the matter and velocity divergence power spectra 
and halo mass function, as well as the halo profiles of the 
scalar degree of freedom, forces, and density. 

5.1 Simulation setup 

All the simulations performed in this study have used the 
same initial conditions, which were generated using 2LPT 
(Crocce et al. 2006) from a ACDM cosmology with Qrn = 
0.269, = 0.731, h = 0.704, ris = 0.966 and erg = 0.8. 

The simulations have N — 512^ particles in a box of size 
B = 250/i“^Mpc and they start at redshift z — 49. All mod¬ 
ified gravity models simulated here have the same expan¬ 
sion history as a ACDM model with the above parameters 
and the evolution of density perturbations at high redshifts 
(z > 10) is almost identical to that of the ACDM model 
justifying the use of the same initial conditions. 

For simulations with a ramses based code, we used a 
coarse-level grid with refinement level /min = 9 correspond¬ 
ing to 512^ coarse cells. Each cell was refined if the num¬ 
ber of particles contained in it exceeded 8. The maximum 
level of refinement corresponded to /max = 15 for F6 and 
TcHo = 5 and /max = 16 for F5 and VcHo = 1. For the 
MG-GADGET simulations the relative tree opening criterion 
was used where the tree is open when the relative force ac¬ 
celeration error is larger than 0.0025, the force softening 
was 18.75 kpc/h and for the long-ranged forces a particle 
mesh grid with 512^ grid cells was used. The force softening 
corresponds to the grid spacing at level ~ 14 in the AMR 
hierarchy. 

The codes compared in this project differ in two aspects: 
(i) the exact way in which the equations of the scalar field 
are solved and (ii) the force calculation and time stepping, 
which were taken from the original codes they were built 
from (e.g. p-gadget3 and ramses). To separate the differ¬ 
ences that arise from these two aspects, simulations of the 
standard ACDM model were also performed. This way, dif¬ 
ferences between the code predictions for the absolute value 
of a measured quantity (e.g. Pk) are affected by both aspects. 
On the other hand, comparisons of the code predictions for 
the relative difference to ACDM (e.g. APfc/T/c,ACDM) should 
be mostly determined by aspect (i) and not so much by as¬ 
pect (ii). 

The run-time of the modified gravity simulations we 
have performed was about 5 — 10 times that of the corre¬ 
sponding ACDM simulation. 

5.2 Matter power spectrum 

As a consistency check, we have measured the matter power 
spectrum, Pk, from the simulations using three independent 
codes. One is the powmes code (Colombi et al. 2009), which 
uses fourier transforms with folding methods to compute Pk . 


Another is a code written by one of the authors of this paper 
(RS), which works similarly to powmes except that it does 
not apply folding methods. Both these codes deconvolve the 
window function of the density assignment and shot-noise is 
subtracted. Lastly, we have also measured the power spec¬ 
trum using a density field obtained with the Delanuay Tesse- 
lation Field Estimator (DTFE) method of Schaap & van de 
Weygaert (2000). We have found that the results from these 
three codes agree very well (to the 1% level ^) and when 
considering the agreement was below 0.1 — 0.5% for 

all scales of interest. All the plots shown in this paper are 
those obtained using the code made by RS mentioned above. 
We show the power-spectrum out to the particle Nyquist 
frequency, /cmax = = 6.4/iMpc“\ of the sim¬ 

ulations. 

Next, we present our results for the matter power spec¬ 
trum, discussing separately the results obtained for ACDM, 
f{R), DGP and Symmetron models. 

5.2.1 ACDM 

Before comparing the results for the modified gravity mod¬ 
els, it is instructive to have a look at how the codes compare 
for ACDM. This is shown in the left panel of Fig. 2 for z = 0 
(solid) and z = 1 (dashed). The red and green lines show, 
respectively, the ratio of the ISIS and egosmog result to 
that of MG-GADGET. One notes that the two RAMSES-based 
codes predict less power (4 — 5% at z = 0,/c ~ 7/iMpc“^) 
than the mg-gadget simulations. This can be attributed 
to the different ways the base codes compute the gravita¬ 
tional force on small scales. In p-gadget 3, the force on large 
scales is computed using a particle-mesh method, just like 
in RAMSES. On small scales, however, p-gadget 3 switches 
to a tree method, whereas in ramses the calculation re¬ 
mains as on large scales. The accuracy of the ramses code 
on small scales depends also on the criteria to trigger a re¬ 
finement of the AMR grid. In all RAMSES-based AMR sim¬ 
ulations of this paper, the grid refines itself whenever the 
particle number inside a given cell exceeds 8. Due to these 
differences between the force calculation on small scales, one 
should therefore not expect perfect agreement between the 
RAMSES-based and GADGET-based codes. We refer the reader 
to Schneider et al. (2015) for a more detailed comparison 
study of the performance of the ramses, gadget and also 
Pkdgrav 3 codes in ACDM simulations. 

The result shown by the blue and pink lines illustrates 
the impact of using different schemes for interpolating the 
density and the force between the particle positions and the 
grid in the two RAMSES-based codes. The blue lines show 
that the ACDM results obtained with isis and egosmog are 
in very good agreement (< 1% error for all scales and times 
shown) if both codes use the same interpolation scheme, in 
this case CIC. This shows that the modifications made to 
RAMSES to develop egosmog and isis do not introduce any 
systematics in the way the codes work for GR On the 
other hand, compared to the egosmog run with TSC, the 

^ For this comparison we have ignored the four largest Fourier 
modes where cosmic variance is significant since the different 
codes use different methods to estimate the power on these scales 
which can lead to quite large differences. 

® We note that the versions of the isis and egosmog used in this 
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power in isis is higher by ~ 3%, (7%) for z = 0, (z = 1), 
for k ^ 7/iMpc“^. This is because the CIC interpolation 
distributes the mass of each particle into fewer grid cells, 
which results in higher peaks in the density field, compared 
to TSC. We note that the size of the differences between 
ISIS and ecosmog to mg-gadget are comparable to the 
differences induced by different interpolation schemes. 


5.2.2 f{R) 

The f{R) model was simulated with egosmog, mg-gadget 
and ISIS. In all codes, the relaxation of the scalar field equa¬ 
tion was performed on an AMR grid. The results are shown 
in the middle and right panels of Fig. 2 for the F5 and F6 
f{R) models, respectively. Fig. 3 shows the relative differ¬ 
ence to ACDM in F5 (left) and F6 (right) for z = 0,1, 2. 
This depicts the known result that the modifications to grav¬ 
ity in this model boost structure formation on small scales 
and that these effects are stronger in the F5 than in the F6 
model. It is perhaps also interesting to note that, for the 
F5 model at z = 0, the relative difference to ACDM does 
not flatten out on scales k < 0.1/iMpc“^, which are scales 
on which the growth is scale-independent in ACDM. This 
shows that in modified gravity, the naive expectation (in¬ 
spired from ACDM) for the scales on which the growth of 
structure is scale-independent can be misleading (Hellwing 
et al. 2013b). 

The differences observed in the right panels of Fig. 2 
are similar in shape and size to those in the left panel for 
ACDM, which suggests that the calculation of the fifth force 
is consistent in between the three codes. This is confirmed 
by the result depicted in Fig. 3, which shows the good agree¬ 
ment between the three codes for all scales and times shown 
In particular, at z = 0 (z = 2), all codes agree to within 
1% for scales k < 7/iMpc“^ (k < 5/iMpc“^). From this, we 
can conclude that any differences between these three mod¬ 
ified gravity codes for f{R) are driven almost exclusively by 
the differences in their main codes (in this case ramses and 
gadget), and not by the extra modules that solve for the 
effects of the fifth force. This is a very reassuring result. 


5.2.3 DGP 

The simulations of the DGP model were performed with 
the dgpm and egosmog codes. The latter was run both 
with refinements and with a fixed grid to better compare 
with dgpm, which is fixed grid only. The top panel of Fig. 4 
shows the ratio of the ACDM power spectrum of the dgpm 
to egosmog simulations without refinements for z = 0 and 
z — 1. The lower panels of Fig. 4 show the relative differ¬ 
ence to ACDM measured in the simulations of the rcHo = 5 
(left) and TcHq = 1 (right) DGP models. We recover the 
known result that, in the DGP model, the amplitude of 


paper are not built on the exact same release of ramses, which 
explains why their ACDM results are not ’exactly’ the same. 

^ Note that mg-gadget agrees well with the two RAMSES-based 
codes on small scales. This shows that any errors arising from 
summation errors in the tree algorithm in gadget do not trans¬ 
late into discrepant power spectrum results. 


the power spectrum is boosted by a scale-independent fac¬ 
tor on scales k < 0.1/iMpc“^. On mildly nonlinear scales, 
0.1/iMpc“^ ^ k < l/iMpc“^, the boost in the power spec¬ 
trum is stronger than on linear scales due to mode-coupling. 
However, on nonlinear scales, k > l/iMpc“^ (halo size 
scales), the suppression effects of the Vainshtein screening 
mechanism are dominant, which effectively reduces the im¬ 
pact of the fifth force on the power spectrum (e.g., Schmidt 
et al. 2010). 

Fig. 4 shows that the three codes agree very well (up 
to 1%) on scales k < l/iMpc“^. For k > l/iMpc“^, how¬ 
ever, the power in the dgpm simulations is higher than in 
egosmog. This is due to the different interpolation schemes 
used in the codes. In particular, in dgpm, the CIC inter¬ 
polation yields a density field with higher peaks than the 
density field in egosmog, which is smoother because of the 
use of the TSC scheme. This is similar to the egosmog and 
ISIS results in the left panel of Fig. 2 (blue and pink lines). 
The lower panels of Fig. 4 also show the egosmog result 
with refinements (blue). For the TcHq = 5 model the three 
codes are, overall, in good agreement for all redshifts and 
scales shown (< 1% for k < 5/iMpc“^). However, the modi¬ 
fications to gravity in the VcHq = 5 model are weaker than 
in the TcHq — 1 case, and as a result, it is easier to inter¬ 
pret the code results for the TcHq — 1 model. For this case, 
the three codes agree very well for k < l/iMpc“^. Note also 
that the egosmog results with refinements agree with its re¬ 
sults for fixed grid on these large scales. For k > l/iMpc“^, 
dgpm is also in good agreement with the results from egos¬ 
mog for fixed grid. Recall that the two codes solve the DGP 
scalar field equation in substantially different ways (cf. Sec¬ 
tion 4.2.2), so this is a nontrivial test. On these small scales, 
the agreement of the fixed grid codes with the egosmog 
code with rehnements gets worse, but this is expected due 
to the gain in resolution in the latter. It is also interesting to 
note that, at z = 0 and for k > 4/iMpc“^, the enhancement 
in the power is smaller in the egosmog simulations with 
refinements compared to the fixed grid cases. The explana¬ 
tion here is that the AMR nature of the grid allows it to 
resolve better the higher density peaks that exist on these 
smaller scales. This results in the code capturing better the 
suppression effects of the screening, and hence, the boost in 
the power relative to ACDM becomes less pronounced. This 
is more pronounced at z = 0 compared to z = 1 because 
at earlier times the density field is less evolved, hence the 
screening efficiency is also weaker. 

In summary, we conclude that the two available N-body 
implementations of the Vainshtein screeening agree well in 
the non-refining case, with the differences from the refined 
case appearing fully consistent with being due to the higher 
resolution of the latter. In the future, it would be desirable 
to also test an independent implementation of Vainshtein 
screening with refinements. 

5 . 2.4 Symmetron 

The left panel of Fig. 5 shows the power spectrum results 
for the Symmetron model, which were obtained with the ISIS 
and ISIS-NONSTATIG code, both run without refinements. The 
figure shows that the impact of the time-derivative terms 
in the equation of the Symmetron model, equation (52), is 
below the 0.5% level for all times and scales shown. More- 
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Figure 2. Comparison of the matter power spectrum results of the ecosmog, mg-gadget and isis codes for the ACDM (left), F5 
(middle) and F6 (right) models, as labelled. When comparing the two RAMSES-based codes (isis and egosmog), we show the egosmog 
results from simulations run with CIC and with TSC interpolation schemes, as labelled. 


over, there seems to be no trend with scale. Hence, we can 
conclude that, in what concerns measurements of the non¬ 
linear matter power spectrum from iV-body simulations of 
the Symmetron model, the use of the quasi-static limit has 
virtually no impact on the results. 

This result is not unexpected since the calculation of 
the matter power spectrum is dominated by high-density re¬ 
gions (haloes, if one thinks about it in the framework of the 
halo model), where the time derivatives are indeed expected 
to be negligible relative to the spatial ones. Consequently, 
it may be of interest to investigate whether the quasi-static 
assumption remains also a good approximation for observ¬ 
ables which are more sensitive to lower density regions. Such 
an investigation is not explored in this study. 

It would also have been good to have a comparison with 
MG-GADGET and EGOSMOG for the symmetron model but 
we leave this to future work. A brief comparison of P(k) 
for the (quasi-static) symmetron model and the codes ISIS, 
EGOSMOG and mlapm can be found in Llinares et al. (2014). 


5.3 Velocity divergence spectra 

We have measured the power spectrum of the velocity di¬ 
vergence field defined as Peeik) = where ^(x) = 

• i’(x), with V being the peculiar velocity field. In the 
linear regime, 0 is related to the matter density contrast as 
0 oc —(5/, where / = dln(5/dlna is the linear growth rate. We 
show only results for the divergence of the velocity field, but 
note that on small scales, where nonlinear processes become 
important, the vorticity (rotational component of r'(x)) is 
nonnegligible and hence the whole velocity field cannot be 
described solely by 0 . 


To measure Pee{k)^ we constructed a volume-weighted 
velocity field (Bernardeau & van de Weygaert 1996) with the 
DTFE method implemented in the publicly available code 
of Cautun & van de Weygaert (2011). We refer the reader to 
Li et al. (2013c) for more details about our method to com¬ 
pute the velocity divergence field from the A-body particle 
positions and velocities. Next, we discuss our results for the 
f{R) and DGP simulations. We have also measured Pee for 
the Symmetron simulations of the ISIS and isis-nonstatig 
codes, but since there are virtually no differences between 
the full and quasi-static results, we refrain from showing 
them. 


5.3.1 f(R) 

Fig. 6 shows the fractional difference of the velocity diver¬ 
gence power spectrum with respect to ACDM in the F5 (left) 
and F6 (right) models. The enhancement in the amplitude 
of Pee in f(R) relative to ACDM is noticeably larger than 
that seen for the matter power spectrum (see also Jennings 
et al. 2012; Li et al. 2013c,b; Hellwing et al. 2014). In partic¬ 
ular, for the F5 (F6) model at z = 0 and k ^ 3/iMpc“^, the 
amplitude of Pee is enhanced by 50% 20%), while the 

boost in the amplitude of P{k) is kept at 20 — 25% (~ 5%) 
only. The velocity field is more sensitive than the density 
field to the modifications to gravity because it starts to be 
affected at earlier times, and the effects get accumulated 
throughout the history of structure growth. 

As in the case of the matter power spectrum, the agree¬ 
ment between egosmog, mg-gadget and isis is notable. In 
particular, for both F5 and F6 at z = 0, the three codes 
agree up to 1% down to /c ~ 3/iMpc“^, and for higher k- 
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Figure 3. Fractional difference of the matter power spectrum with respect to ACDM from the simulations of the F5 (left) and 
F6 (right) models performed with the ecosmog (TSC), mg-gadget and isis codes, as labelled. In the lower panel insets, e = 
(P/Pacdm) code / (^/^ACDM)ref “ with MG-GADGET being the reference code. At 2 : = 0, the codes are all accurate to within 1% 
on all scales shown. 


values the difference never exceeds 2%. The case of the ISIS 
code for F5 at z = 1 is an exception to this very good agree¬ 
ment, for which the difference with the other two codes is at 
the 2 — 4% level for k ~ 0.3 — 3/iMpc“^. We have checked 
that at 2 : = 2 the agreement between the three codes is at 
the 1 — 2% level, which suggests that the larger disagree¬ 
ment at z = 1 could be due to a transient effect related, for 
instance, to differences in the time integration We leave 
a more thorough investigation of this difference for future 
work. 


5.3.2 DGP 

Fig. 7 shows the fractional difference of the velocity di¬ 
vergence power spectrum with respect to ACDM in the 
TcHq — 5 (left) and TcHq — 1 (right) DGP models. The 
agreement between dgpm and ecosmog is very good, with 


® For example, slight differences in the AMR structure and ve¬ 
locity field in the egosmog and isis simulations can result in 
different time step sizes (determined by the same criteria as in 
standard ramses). 


any differences in TcHq — 5 predictions being below ~ 1% 
for all times and scales shown. For the rcHo = 1 model the 
agreement between the codes worsens, but the difference is 
always below the 2% level. Similarly to the case of the f{R) 
model, we also note that the modifications to gravity in the 
DGP model affect the amplitude of Pee more than they af¬ 
fect the amplitude of P{k). 

In Fig. 7, it is interesting to note that on small scales 
{k > 2 — 3/iMpc“^), the difference between the fixed (red 
and green) and refined grid (blue) results is smaller than 
that seen in Fig. 4 for the matter power spectrum. Given 
the gain in resolution when ecosmog is run with refine¬ 
ments, one does not expect the results to fully agree with 
fixed grid simulations which cannot resolve small scale struc¬ 
tures, and as a result, the agreement depicted in Fig. 7 may 
seem surprising. We do not perform any detailed investiga¬ 
tions of this result, but simply note that the density field 
used to compute P(k) is mass-weighted, whereas the den¬ 
sity field used to compute Pee is volume-weighted. This, to¬ 
gether with the suppressed contribution of virial velocities 
to the velocity divergence, may help explain why the use of 
adaptively refined grids does not have a critical impact on 
the resulting Pee- For completeness, we further note that. 
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Figure 4. The top panel shows the comparison of the matter power spectrum results of the dgpm and ecosmog (fixed grid) simulations 
of the ACDM, TcHq = 5 and TcHo = 1 DGP models, as labelled. The lower left and lower right panels show the fractional difference 
with respect to ACDM of the two codes for the VcHo = 5 and TcHo = 1 DGP models, respectively. The lower panels show the egosmog 
results both for fixed and refined grid simulations, as labelled. In the lower panel insets, e = {P/Pacom) code / (P/Pacom)— 1, with 
EGOSMOG (refined grid) being the reference code. 


as found in Falck et al. (2015), peculiar velocities on small 
scales are also less affected by the Vainshtein mechanism, 
compared to the effects of the Chameleon mechanism. 


5.4 Halo mass function 

The halo mass function, n(M), defined as the number den¬ 
sity of dark matter haloes with mass M, is an important 
statistic that is particularly sensitive to modifications to 
gravity. We used a modified version of the spherical over- 
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Figure 5. Comparison of the matter power spectrum (left) and halo mass function (right) results of the Symmetron model simulations 
performed by the isis and isis- nonstatic codes. The lower panels shows the error induced by applying the quasi-static approximation, 
uasi-static/^Full “ I5 whcreas the upper panels shows the fractional difference relative to ACDM, as labelled. In all panels, the 
error induced by employing the quasi-static limit lies comfortably below 1% for all the times and scales shown. 


density Amiga’s Halo Finder (ahf) code (Gill et al. 2004; 
Knollmann & Knebe 2009) to identify dark matter haloes 
and calculate their profiles. For this, simulation outputs from 
all codes were written in the standard gadget2 format, but 
with new data blocks added, which contain the components 
of standard gravity and the fifth force as well as the New¬ 
tonian potential and the scalar field at particle positions. 
We outputted quantities at particle positions rather than 
leaf cells of the AMR grids or trees because, unlike the leaf 
cells, the particle IDs are the same in all simulations, mak¬ 
ing comparisons more straightforward. Note that this means 
that the scalar field and forces in underdense regions were 
not very well sampled, which however should not be a seri¬ 
ous issue given that we are mostly interested in haloes. 

The major modifications to ahf were threefold: (i) new 
routines to read the above data, which has the same format 
as the default gadget2 data such as particle coordinates 
and velocities; (ii) new routines to compute the force and 
scalar field profiles in haloes, by averaging over their values 
at all particle positions in a given spherical shell (the radial 
binning scheme for this was the same as in the default ahf 
code); (iii) the routine in ahf which does the removal of 
unbound particles was also modified so that the code deter¬ 
mined whether a particle was bound or not by comparing its 
velocity with the total gravitational potential instead of the 


standard Newtonian potential. Note that if (iii) is not prop¬ 
erly done, then the mass function tends to be lower because 
more particles are considered as too fast to be bound. Li 
& Zhao (2010) studied the effect of taking account the fifth 
force in the halo unbinding process for certain chameleon 
models using A-body simulations and found a noticeable 
difference in the resulting mass functions (see also Hellwing 
et al. (2010, 2013a)). 

We used ahf with Avir = 200 so that Mahf = M2ooc- 
Next, we discuss our mass function results for the f{R), 
DGP and Symmetron simulations. 

5.11 f{R) 

Figs. 8-9 show the mass function results of the egosmog, 
ISIS and mg-gadget simulations for f{R). As first shown 
quantitatively in Schmidt et al. (2009), f{R) models pre¬ 
dict an enhancement in the abundance of haloes relative 
to AGDM. The enhancement is more pronounced in the F5 
model because the Ghameleon screening is less efficient. This 
is particularly noticeable at the high-mass end for which, in 
the F6 model, the number density of haloes is almost the 
same as in AGDM. The differences in the mass dependence 
of An/riACDM for different f{R) model parameters and red- 
shifts illustrates the complex interplay between the mass 
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Figure 6. Fractional difference of the velocity divergence power spectrum with respect to ACDM from the simulations of the F5 
(left) and F6 (right) models performed with the ecosmog (CIC), mg-gadget and isis codes, as labelled. In the lower panel insets, 
e = (F/FACDM)code / (^/^ACDM)ref “ with MG-GADGET being the reference code. 


and time dependence of the Chameleon screening mecha¬ 
nism (see e.g. Lombriser et al. 2014; Shi et al. 2015; Gronke 
et al. 2015a; Gronke et al. 2015b, 2014; Winther et al. 2012, 
for studies of halo properties in f{R)). 

Fig. 8 shows the ratio of the ecosmog and isis results to 
those of MG-GADGET for the ACDM (left panel), F5 (middle 
panel) and F6 (right panel) models. The three codes show 
varying levels of agreement (between 2-10%) throughout the 
mass range probed by our simulations and for the two red- 
shifts shown. In particular, at 2 : = 0 and intermediate mass 
scales, M ~ 10^^ MqIH, ecosmog and isis agree with mg- 
gadget at < 4%, but the agreement worsens to 5 — 10% for 
smaller mass scales. At the high-mass end, M ~ 10^^ M^jh, 
the finite box size limits us to only a few halo samples, which 
is why the high mass end is noisy (as can be checked by the 
scatter between different mass bins). Overall, the trend is for 
the two RAMSES-based codes to underpredict the abundance 
of smaller mass haloes, compared to mg-gadget. This dis¬ 
crepancy can be linked to the way the base codes ramses 
and p-gadget3 compute the force on small scales (see dis¬ 
cussion in Sec. 5.2). Also, there is a marked improvement 
in the agreement between ecosmog and ISIS when the two 
codes are run with the same density assignment scheme (GIG 
in this case). 

The amplitude and shape of the curves in Fig. 8 is simi¬ 
lar for AGDM and the two f(R) models, which indicates that 


differences between the modified gravity codes are mostly 
driven by differences in the default codes. Indeed, this is 
again confirmed by the result of Fig. 9 which shows the 
fractional difference in the f(R) mass functions relative to 
AGDM. One notes that the three codes agree very well, es¬ 
pecially at 2 ; = 0, echoing the results seen in the previous 
two sections for the matter and velocity divergence power 
spectra. 

5 . 4.2 DGP 

The top panel of Fig. 10 shows the ratio of the mass func¬ 
tion results obtained with the dgpm code to those obtained 
with ECOSMOG (both with and without refinements) for the 
simulations of the AGDM and VcHo = 5 and VcHo = 1 DGP 
models. For the three models, there is a systematic trend for 
DGPM to produce more massive haloes than ecosmog, even 
when ECOSMOG is run without refinements like dgpm. This 
result can be linked to the different density assignments of 
the two codes (GIG for dgpm vs TSG for ecosmog). The 
lower panels of Fig. 10 show that the agreement between the 
two codes improves when one looks at the fractional differ¬ 
ence of the two DGP models to AGDM. For the case of the 
TcHq = 5 model, although there are still visible differences 
between the results of the two codes, these remain of the 
same order as the bin-to-bin scatter. The agreement wors- 
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Figure 7. Fractional difference of the velocity divergence power spectrum with respect to ACDM from the simulations of the VcHo = 5 
(left) and VcHo = 1 (right) DGP models performed with the dgpm and ecosmog codes. The two sets of egosmog results correspond to the 
results from simulations run with fixed and refined grids, as labelled. In the lower panel insets, e = (P/Pacdm) code / (^/^ACDM)ref “ ^5 
with EGOSMOG (refined grid) being the reference code. 


ens slightly for the rcHo = 1 model. Overall, both dgpm 
and EGOSMOG are in good agreement in their predictions 
for An/nACDM, although to a lesser extent than the agree¬ 
ment for the matter and velocity power spectra seen in the 
previous sections (Figs. 4 and 7). 

5.4-3 Symmetron 

The right panel of Fig. 5 measures the impact of the quasi¬ 
static limit on the mass function of the Symmetron model, 
as predicted by the ISIS and isis-nonstatig codes (both run 
without rehnements). As seen in the case of the mass power 
spectrum (upper panels), relaxing the quasi-static approxi¬ 
mation has no appreciable effect on halo abundances: differ¬ 
ences are < 0.5% and have no clear dependence on mass. 

5.5 Halo profiles 

Finally, we turn our attention to the code results for the 
radial prohles of the density, velocity dispersion, force and 
scalar held around dark matter haloes. For this, we used 
the z = 0 AHF halo catalogues of each simulation and 
binned them in mass according to: [5 x 10^^, 1 x 10^^], 
[1 X 10^^ 5 X 10^^], [5 X 10^^ 1 X 10^^] and [1 x 10^^, 5 x 10^^] 
MqIH. For each mass bin, we scaled the ahf haloes by their 


virial radii (determined by the ahf code) and stacked them 
to compute average prohles of the different quantities. We 
take the variance of this average as the errorbars. The den¬ 
sity, scalar held, and force prohles are calculated using ah 
particles around a given ahf halo centre, which includes par¬ 
ticles that he beyond the ahf halo. On the other hand, the 
velocity dispersion prohles are calculated using only parti¬ 
cles dehned to be part of the ahf halo, thus they only extend 
out to the virial radius. 

As before, we discuss our results for the halo prohles in 
turn for the f{R), DGP and Symmetron models. 


5.5.1 f(R) 

We plot the prohles of the scalar held in Fig. 11, of the force 
modulus in Fig. 12, and of the density and velocity disper¬ 
sion in Fig. 13 for the haloes found in the F5 and F6 egos¬ 
mog, MG-GADGET and ISIS simulations. The egosmog re¬ 
sults in these hgures correspond to the runs performed with 
the CIC interpolation scheme (as in isis). For mg-gadget, 
the force prohles shown are those calculated by interpolating 
the gradient of the scalar held from the grid to the particle 
positions instead of the effective density method (recall the 
discussion about Fig. 1 in 4.2.1). The gradient method cap- 
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Figure 8. Comparison of the halo mass function results of the ecosmog, mg-gadget and isis codes for the ACDM (left), F5 (middle) 
and F6 (right) models, as labelled. The two sets of egosmog results correspond to the results from simulations using the CIC and TSC 
density assignment, as labelled. 


tures more accurately the suppression effects of the screening 
in the inner regions of the haloes. 

Overall, the ISIS and ecosmog codes agree very well, 
with any deviations typically lying within the error bars in 
the profiles depicted in Figs. 11, 12 and 13. This is reassuring 
but not very surprising, considering that they are both based 
on RAMSES and are run with the same settings. It is there¬ 
fore more interesting to compare ISIS and ecosmog with 
MG-GADGET, for which some differences exist. For instance, 
although at large radii there is good agreement between the 
scalar field profiles obained by the three codes, in the inner 
regions of the haloes there is not. This is particularly notice¬ 
able in the F6 model at r/rvir ^ 1 (right panels of Fig. 11), 
for which MG-GADGET ovcrpredicts the values of fn com¬ 
pared to ECOSMOG and isis. In the F5 model the discrepan¬ 
cies soften considerably, especially for the low mass haloes. 
These differences in the fn profiles naturally translate into 
differences in the amplitude of the fifth force on small radial 
scales, as seen in the lower four panels of Fig. 12. However, 
the lower panels of Fig. 12 show that the discrepancies be¬ 
tween MG-GADGET and the two RAMSES-based codes are only 
appreciable when the fifth force is a small fraction (< 0.1) 
of the total force (due to screening). This means that, al¬ 
though the codes may disagree on their exact predictions 
for the amplitude of the modifications to gravity, they only 
do so in regimes where the fifth force is not very important 
anyway. This is confirmed by the fact that the density and 
velocity dispersion profiles shown in Fig. 13 agree very well 
(there are differences of order 5%, which we note are likely 
to come from differences in the base codes.) 

At large radii (r/rvir ^ 1), mg-gadget overpredicts the 
amplitude of the Newtonian force compared to egosmog 
and ISIS. Here, we stress that what we average in a given 
halo mass bin is the force modulus and not any directional 
component of the force (e.g. radial). This is why the force 


profiles do not keep decaying torwards large radii but instead 
level off due to the matter distribution that surrounds the 
haloes^. This can explain the mismatch in the Newtonian 
force at large radii because the matter distribution around 
haloes in gadget and ramses simulations is not exactly 
the same, and the base code algorithms that compute the 
Newtonian force are also different. What is important here 
is that, at these large radii, the three codes agree very well 
in their fifth force predictions, which are the corrections to 
normal gravity we are interested in testing in this paper. 

5.5.2 DGP 

Figs. 14 and 15 show the same as Figs. 12 and 13, but for the 
DGP simulations performed with the dgpm and egosmog 
(with and without refinements) codes. The absolute value 
of the scalar field in the DGP model is irrelevant as the 
equations of the model contain only its derivatives. For this 
reason, we do not show the scalar field profiles and prefer to 
plot the gradient of the field (which, up to a factor of 1/2, 
is the fifth force as seen in equation 15). 

The Newtonian force profiles of the two codes are in very 
good agreement when egosmog is run without refinements 
(like dgpm). There are marked differences in the solutions of 
the two codes for the size of the fifth force in the inner regions 
of the haloes, but since the amplitude of the fifth force is 
small there anyway, the difference does not translate into the 
density and velocity dispersion profiles (as seen in Fig. 13). 
However, when analysing these results, it is important to 
bear in mind that the grid size of the nonrefined simulations 

® If one would average the radial component of the force for a 
large number of haloes, then the contribution from the surround¬ 
ing structure would cancel out. 
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Figure 9. Fractional difference of the halo mass function with respect to ACDM from the simulations of the F5 (left) and 
F6 (right) models performed with the ecosmog (TSC), mg-gadget and isis codes, as labelled. In the lower panel insets, e = 
(n/nACDM)code / (^/'^ACDM)ref “ With MG-GADGET being the reference code. 


(Ar ~ 0.5/i“^Mpc) is half of the typical halo size scales of 
l/i“^Mpc. As a result, one should not attempt to draw any 
physically meaningful conclusions from the results depicted 
by the red and green lines. We showed these lines simply to 
illustrate that the dgpm and ecosmog (nonrefined) codes 
agree in their total force profiles, even though their solutions 
may not be accurate. 


When ECOSMOG is run with refinements on the grid, 
the code is able to better resolve the matter distribution 
in and around the haloes. The top panels of Fig. 15 show 
that the density profiles of the nonrefined simulations only 
agree with those of the refined simulation for r/rvir ^ 5 — 10. 
This provides a measure for the radial scales below which 
one should not trust the nonrefined simulations. It is also 
interesting to note that the ratio of the fifth to Newtonian 
force does not depend on halo mass, as seen in the lower 
panels of Fig. 14. This illustrates that the efficiency of the 
Vainshtein mechanism in the DGP model is independent 
of the mass of the haloes (Schmidt 2010; Falck et al. 2014, 
2015), which is different from what is seen in the lower panels 
of Fig. 12 for the f{R) models. 


5.5 .3 Symmetron 

Fig. 16 shows the radial profiles of the scalar field and of the 
fifth to Newtonian force ratio obtained from the ISIS and 
isis-NONSTATic simulations of the Symmetron model. The 
scalar field in the isis-nonstatic simulations oscillates very 
rapidly with time. For this case, the profiles we show repre¬ 
sent the mean value averaged over several oscillations close 
to z = 0. We note also that since these simulations were per¬ 
formed on a fixed grid, they suffer from the same resolution 
problems as the non-refined DGP runs. This prevents the 
simulations from fully capturing a number of effects such 
as the efficiency of the screening. Nevertheless, for the sake 
of determining the impact of the quasi-static limit these is¬ 
sues can be ignored. The result of Fig. 16 reinforces the 
conclusions drawn previously that relaxing the quasi-static 
limit has little impact on the simulations of the Symmetron 
model. 


6 SUMMARY AND CONCLUSIONS 

A-body simulations of modified gravity play a key role in 
cosmological tests of gravity using large-scale structure ob- 
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Figure 10. The top panel shows the comparison of the halo mass function results of the dgpm and ecosmog (fixed grid) simulations 
of the ACDM, TcHq = 5 and TcHq = 1 DGP models, as labelled. The lower left and lower right panels show the fractional difference 
with respect to ACDM of the two codes for the TcHq = 5 and TcHq = 1 DGP models, respectivly. The lower panels show the egosmog 
results both for fixed and refined grid simulations, as labelled. In the lower panel insets, e = (?T'/'^^'ACDM)code / (^/'^ACDM)ref “ ^5 with 
EGOSMOG (refined grid) being the reference code. 


servations. In this work, we have performed an extensive 
comparison of the simulation results of /(i?), DGP and Sym- 
metron gravity using five modified gravity iV-body codes, 
which were the dgpm, egosmog, mg-gadget, isis and isis- 
NONSTATIG codcs. In the models we simulated, the gravita¬ 


tional law is modified due to the presence of a fifth force 
mediated by a scalar field. The algorithms in these codes 
differ from those for standard gravity by having extra mod¬ 
ules that iteratively relax the equation of the scalar field on a 
grid/mesh to determine the fifth force. The isis-nonstatig 
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Figure 11. Radial profiles of the fji scalar field in dark matter haloes found in the simulations of the ecosmog (CIC), mg-gadget and 
ISIS codes for the F5 (left) and F6 (right) models. The result is shown for the following four halo mass bins: M G [1 x 5 x 10^"^], 

M G [5 X 10^^, 1 X 10^"^], M G [1 X 10^^, 5 X 10^^] and M G [5 x 10^^, 1 x 10^^] Mq/Zi (from bottom to top). For clarity, the profiles for 
each mass bin are displaced vertically, with the four horizontal lines (also displaced vertically) indicating Jr/Zro = 1 of each mass bin. 


code is a version of the ISIS code that goes beyond the 
quasi-static approximation that is often employed in mod¬ 
ified gravity studies. We used this code to test the validity 
of this assumption in the case of the Symmetron model. 

The modified gravity routines included in the mod¬ 
ified gravity codes are typically installed in existing and 
well tested iV-body codes for GR (such as ramses or p- 
gadgetS). Since these main codes can show some level 
of disagreement in their results for standard ACDM, then 
these differences would naturally propagate into the modi¬ 
fied gravity results. To overcome this, in this paper we have 
run also simulations of standard ACDM, which were used 
as a reference to measure the effects of the modifications to 
gravity. Hence, our main concern in this paper was to com¬ 
pare the code predictions for the fractional difference with 
respect to ACDM of a given quantity, and not so much its 
absolute values. 

All our simulations of the different models and codes 
start from the same set of initial conditions, which allows 
a direct comparison between the results. We analysed the 
code results for the power spectrum of the matter density 
fluctuations and peculiar velocity divergence. We have also 
compared results for the abundances of dark matter haloes 
and for their density, velocity dispersion and force profiles. 

We have generally found agreement at the few-pereent 
level in the properties of the matter density and veloeity 
fields. This means that modified gravity simulations satisfy 
the aecuracy requirements of eurrently planned large-seale 
strueture surveys, provided that the absolute calibration of 
GR predictions is of comparable accuracy. 

In what follows, we recap the main results of this com¬ 
parison project in more detail. 


• Matter power spectrum 

Given its immediate relation to galaxy redshift and lensing 
observables, the matter power spectrum P{k) is one of 
the most important observables for tests of gravity. For 
AGDM, we found that the modified gravity codes (run 
with the routines for modified gravity switched off) agree 
up to 1% for k < l/iMpc“^, but start differing on smaller 
scales (~ 5% at /c ~ 5/iMpc“^) due to different density 
and force assignments schemes and intrinsic algorithmic 
differences in the force calculation in AMR (for egosmog 
and isis) and TreePM (for mg-gadget) codes (cf. Fig. 2). 
For the f{R) simulations, which were performed with the 
EGOSMOG, MG-GADGET and ISIS codes, we find that any 
differences in P(k) are driven almost exclusively by the 
differences in the base codes (ramses and p-gadget3 
respectively). In terms of the relative difference to AGDM, 
all code results for f{R) agree to better than 1% for k < 
7/iMpc“^ (cf. Fig. 3). 

While EGOSMOG and isis are based on the same GR 
code and use very similar algorithms, mg-gadget is suffi¬ 
ciently different to make this a nontrivial consistency test. 
The simulations of the DGP model were performed with 
the DGPM and egosmog codes. The former code does not 
have an adaptive mesh, but egosmog can be run on both 
a fixed and refined mesh. For the fixed grid case, the two 
codes are in very good agreement in their predictions for 
the relative difference to AGDM (< 1% for /c < 7/iMpc“^) 
(cf. Fig. 4). This is a nontrivial test since the DGP scalar 
field solvers in dgpm and egosmog are substantially dif¬ 
ferent. The fixed grid simulations start to differ from the 
refined grid simulations on small scales, but this is ex- 
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pected since the latter is able to resolve small scales much 
better. 

For the test of the validity of the quasi-static limit in 
the Symmetron model, we have seen that the results 
from ISIS and isis-nonstatic are nearly indistinguishable 
(cf. Fig. 5). We conclude that the impact of adding the 
time derivative terms to the Symmetron field equation of 
motion is negligible (< 0.5% on all scales and redshifts 
shown). We note, however, that in the current implemen¬ 
tation of the time derivatives in the isis-nonstatic code, 
the particles do not feel the details of the rapid oscillations 
of the scalar field, due to the difference in the particle and 
scalar field time steps (see Llinares & Mota 2014, for more 
details). In the future, it would be of interest to clarify 
whether using the same time steps for particles and the 
scalar field has an impact on the results presented here. 

• Velocity power spectrum 

In terms of the velocity divergence power spectrum, Pee^ 
we have found that ecosmog, mg-gadget and isis are 
again in very good agreement in their predictions for the 
fractional difference with respect to ACDM in f{R) mod¬ 
els. At z = 0, the three codes agree to better than 1% for 
k < 3/iMpc“^, for both the F5 and F6 models (cf. Fig. 6). 
For the DGP simulations, the dgpm and ecosmog codes 
are also in very good agreeent, with their predictions dif¬ 
fering by < 2% for both models simulated and for all times 
and scales shown (cf. Fig. 7). Interestingly, we have also 
seen that the Pee results of the DGP simulations on small 
scales are not critically affected by the use of a fixed or 
refined grid. For brevity, we did not show the impact of 
the quasi-static limit on Pee in the Symmetron model, but 
we have checked that it is negligible (just like for P{k)). 
These results are encouraging as the precise modelling of 
cosmic velocity fields is a crucial ingredient in connect¬ 
ing the theoretically predicted clustering statistics in real 
space with the observed galaxy redshift space power spec¬ 
trum. 

• Halo mass function 

The good agreement outlined above for the matter and pe¬ 
culiar velocity divergence power spectrum holds also for 
the halo mass function. Given the relatively small sim¬ 
ulation volume, the accuracy of the mass function com¬ 
parison is limited by cosmic variance. However, we are 
still able to conclude that the codes agree to a satisfac¬ 
tory level. In particular, for f{R), the code predictions for 
An/nACDM agree up to 4% for all mass scales probed 
by our simulation box (cf. Fig. 9). For the DGP simula¬ 
tions the agreement is at the 5 — 10% level (cf. Fig. 10). 
Like for P(k), the inclusion of the time derivative terms 
in the Symmetron model equations leads to negligible dif¬ 
ferences only (< 0.5%). 

• Halo profiles 

Understanding dark matter halo profiles in modified grav¬ 
ity is crucial to devise tests of gravity using galaxy clus¬ 
ters (Lam et al. 2012; Terukina & Yamamoto 2012; Lam 
et al. 2013; Zu et al. 2014; Lombriser et al. 2012; Wilcox 
et al. 2015; Barreira et al. 2015; Terukina et al. 2015). 
For the f{R) simulations, although there are significant 
differences in the fifth force profiles of the RAMSES-based 


codes and mg-gadget in the inner regions of the haloes 
(cf. Fig. 12), these do not translate into differences in the 
resulting density and velocity divergence profiles because 
they only appear in regions where the force modification 
is highly suppressed by chameleon screening (cf. Fig. 13). 
Similarly, the fifth force halo profiles of the dgpm and 
EGOSMOG (fixed grid) DGP simulations also differ, but 
only in the inner regions of the haloes, where the fifth 
force is already weak (cf. Fig. 14). We note, however, that 
although it is reassuring that the two codes agree in their 
halo density and velocity dispersion profiles (which is a 
cross check of the validity of the algorithms), the resolu¬ 
tion of the fixed grid prevents us from trusting their phys¬ 
ical results on halo size scales. Indeed, the DGP simula¬ 
tions run with EGOSMOG (refined grid) show very different 
results, due to the gain in resolution (cf. Fig. 15). 

The halo profiles in the Symmetron simulations performed 
with the ISIS and isis-nonstatig codes are nearly indis¬ 
tinguishable, which reinforces further the validity of the 
quasi-static limit (cf. Fig. 16). 

Tests of gravity on large scales using galaxy power spec¬ 
tra, galaxy dynamics, cluster abundance, cluster profiles and 
so on are one of the main science drivers of many upcoming 
large surveys. Accurate simulations of nonlinear structure 
formation in these models are therefore crucial to ensure ro¬ 
bust theoretical predictions that can be compared with the 
observational data. In this paper, we have seen that the N- 
body codes that have been developed so far are in very good 
agreement for different models. This dedicated comparison 
project constitutes an important validity check of the dif¬ 
ferent codes, which brings us one step closer to performing 
larger and more expensive modified gravity simulations to 
be used to prepare for several future observational efforts. 
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Figure 12. Radial profiles of the Newtonian force (upper), fifth force (middle) and fifth to Newtonian force ratio (bottom) for haloes in 
the F5 (left) and F6 (right) simulations, performed with the ecosmog (CIC), mg-gadget and isis codes, as labelled. The result is shown 
by splitting haloes into the same mass bins used in Fig. 11. In this case, however, the results for each bin are not displaced vertically. 
In the lower panels, the horizontal lines show Ff^/F^ = 1/3, which is the expected value in the absence of screening. When determining 
the force profiles in each mass bin, what is averaged is the force modulii of the haloes, and not its radial component. 
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Figure 13. Radial profiles of the density (upper) and velocity dispersion (bottom) for haloes in the F5 (left) and F6 (right) simulations, 
performed with the ecosmog (CIC), mg-gadget and isis codes, as labelled. The result is shown by splitting the haloes into the same 
mass bins used in Figs. 11 and 12. Note that in the lower panels the profile is only shown up to r/rvir = 1- 
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Figure 14. Radial profiles of the Newtonian force (upper), fifth force (middle) and fifth to Newtonian force ratio (bottom) for haloes in 
the TcHq = 5 (left) and TcHq = 1 (right) simulations, performed with the dgpm and ecosmog codes. The egosmog results are shown 
for both the fixed and refined grid simulations, as labelled. The result is shown by splitting haloes into the same mass bins as used in 
Fig. 12. When determining the force profiles, what is averaged is the force modulii, and not its radial component. Note also that for all 
mass scales, the AMR nature of the grid plays a key role in the measured force profiles. 
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Figure 15. Radial profiles of the density (upper) and velocity dispersion (bottom) for haloes in the VcHo = 5 (left) and rcHo = 1 
(right) simulations, performed with the dgpm and ecosmog codes. The egosmog results are shown for both the fixed and refined grid 
simulations, as labelled. The result is shown by splitting the haloes into the same mass bins used in Fig. 14. For clarity, the result for 
the different mass bins has been displaced vertically in the density panels. In these, the differences between the two codes for fixed grid 
are kept below 5% for all mass bins and scales shown. Note that in the lower panels the profile is only shown up to r/rvir = 1- 
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Figure 16. Radial profiles of the scalar field (left) and of the fifth to Newtonian force ratio (right) for haloes in the Symmetron model 
simulations, performed with the isis and isis-nonstatic codes. The result is shown by splitting the haloes into the same mass bins used 
in the figures of the f{R) and DGP models. As in Fig. 5, the quasi-static limit remains an extremely good approximation. 
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